Modern Analogue Technique

WA take a species approach to reconstruction — each species in the fossil sample that is also in the training set contributes to the reconstructed values

MAT takes a more holistic approach — we predict on basis of similar assemblages

In MAT, only the most similar assemblages contribute to the fitted values

MAT is steeped in the tradition of uniformitarianism — the present is the key to the past

We take as our prediction of the environment of the past, the (possibly weighted) average of the environment of the \(k\) sites with the most similar assemblages

Several things to define; \(k\), (dis)similarity

MAT is \(k\) nearest neighbours (\(k\)-NN) regression/calibration

Modern Analogue Technique

If you want to fit MAT models, my J. Stat. Soft. paper has fully-working examples to follow

Simpson, G.L., 2007. Analogue Methods in Palaeoecology: Using the analogue Package. J. Stat. Softw. 22, 1–29.

Or use rioja::MAT(), especially if you want h-block CV

rioja

analogue is getting a little out-dated now — it needs some love

rioja by is very similar, provides more functionality (WA-PLS, Maximum likelihood TFs, …) but most importantly, h-block cross-validation!

analogue has better support for the Modern Analogue Technique

library("rioja")
wa_r <- WA(y = ik, x = SumSST)
wa_r

Method : Weighted Averaging
Call   : WA(y = ik, x = SumSST) 

Tolerance DW       : No 
Monotonic deshrink : No 
No. samples        : 61 
No. species        : 27 
Cross val.         : none 

Deshrinking regression coefficients:
      Inverse d/s Classical d/s
wa.b0     -5.6876        5.8905
wa.b1      1.2659        0.7246

Performance:
          RMSE      R2  Avg.Bias  Max.Bias    Skill
WA.inv  2.0188  0.9173         0    3.8155  91.7299
WA.cla  2.1078  0.9173         0    3.1324  90.9842

WA-PLS

Principal Component Regression — decompose \(\mathbf{Y}\) into PCA axes and use the first \(l\) of these axes to predict \(\mathbf{x}\): analogue::pcr()

Partial Least Squares is similar — find orthogonal components in \(\mathbf{Y}\) that are most correlated with \(\mathbf{x}\), and use the first \(l\) to predict \(\mathbf{x}\)

WA-PLS is a non-linear version of PLS — use row & column sums to turn a linear method into a non-linear one, just like CCA does: rioja::WAPLS()

An important aspect of fitting these models is to choose \(l\) — the number of components to retain

Bias

Bias is the tendency for the model to over or under predict

Average bias is the mean of the residuals

Maximum bias is found by breaking the range of the measured environment into \(n\) contiguous chunks (\(n = 10\) usually)

Within each chunk calculate the mean of the residuals for that chunk

Take the maximum value of these as the maximum bias statistic

Bias

op <- par(mfrow = c(1,2))
plot(wa_m)

par(op)

Crossvalidation

Without cross-validation, prediction errors, measured by RMSEP, will be biased, often badly so, because we use the same data to both fit & test the model

Ideally we’d have such a large training set that we can split this into a slightly smaller training set and a small test set

Palaeoecological data is expensive to obtain — in money and person-hours!

Also these ecosystems are complex, species rich, noisy etc., so we want to use all our data to produce a model

One solution to this problem is to use cross-validation

General idea: perturb the training set in some way, build a new model on the perturbed training set and assess how well it performs

If we repeat the perturbation several time we get an idea of the error in the model

Several techniques; n-fold, leave-one-out (LOO), bootstrapping, h-block CV

Crossvalidation in analogue

In analogue, several methods are available

For MAT models, LOO is built into the procedure so only bootstrapping is available

For WA models, both LOO and bootstrapping currently available

\(n\)-fold CV & h-block will be available in a future version (!)

Leave-one-out CV

LOO CV is very simple

In turn, leave out each sample from the training set

Build a model on the remaining samples

Predict for the left out sample

Calculate the RMSEP of these predictions

Leave-one-out CV

loo_pred <- predict(wa_m, V12.122, CV = "LOO", verbose = TRUE)

Leave one out sample 10 
Leave one out sample 20 
Leave one out sample 30 
Leave one out sample 40 
Leave one out sample 50 
Leave one out sample 60 

analogue::performance(wa_m)

    RMSE       R2 Avg.Bias Max.Bias 
   2.019    0.917    0.000   -3.815 

analogue::performance(loo_pred)

   RMSEP       R2 Avg.Bias Max.Bias 
   2.218    0.900   -0.014   -4.599 

Leave-one-out CV

rioja::crossval(wa_r, verbose = FALSE)

Method : Weighted Averaging
Call   : WA(y = ik, x = SumSST) 

Tolerance DW       : No 
Monotonic deshrink : No 
No. samples        : 61 
No. species        : 27 
Cross val.         : loo 

Deshrinking regression coefficients:
      Inverse d/s Classical d/s
wa.b0     -5.6876        5.8905
wa.b1      1.2659        0.7246

Performance:
               RMSE      R2  Avg.Bias  Max.Bias    Skill
WA.inv       2.0188  0.9173    0.0000    3.8155  91.7299
WA.cla       2.1078  0.9173    0.0000    3.1324  90.9842
WA.inv_XVal  2.2179  0.9003   -0.0137    4.5985  90.0176
WA.cla_XVal  2.3073  0.9018   -0.0067    3.9326  89.1969

n-fold CV

n-fold CV or leave-group-out CV

Same as LOO but instead of leaving out one sample at a time, we leave out an entire group of samples

n is usually 10

h-block CV

h-block CV very similar to LOO & \(n\)-fold

Instead of breaking training set into n groups & leaving one out at a time, we leave out all observations withinh distance of the target sample

Repeat for each sample in turn (like LOO)

Useful when training data are autocorrelated — as they often are in marine settings

USe rioja::crossval() to do h-block

Bootstrapping I

Bootstrapping used in machine learning to improve predictions

Use bootstrapping to get more realistic RMSEP and bias statistics

We draw a bootstrap sample (sampling with replacement) of the same size as our training set

Build a model on the bootstrap samples

Predict for the out-of-bag (OOB) samples

Bootstrapping II

Bootstrap prediction for each model sample is the mean of the OOB prediction for each sample

Calculate the residuals and then the RMSEP

\[\mathrm{RMSEP_{boot}} = \sqrt{s_1^2 + s_2^2}\]

\(s_1^2\) is the standard deviation of the OOB residuals

\(s_2^2\) is the mean of the OOB residuals

We can also calculate the more usual RMSEP \(\sqrt{\sum_{i=1}^n (y_i - \hat{y}_i)^2 / n}\)

Bootstrapping II

set.seed(1234)
sst_boot <- bootstrap(wa_m, n.boot = 999, verbose = FALSE)
sst_boot

    Weighted Averaging Bootstrap Predictions

Call:
bootstrap(object = wa_m, n.boot = 999, verbose = FALSE) 

Deshrinking     : Inverse 
Crossvalidation : bootstrap 
Tolerance DW    : No 

Performance:
   RMSEP        R2  Avg.Bias  Max.Bias  
  2.2919    0.9012   -0.0298   -4.6312  

Training set predictions:
  V14.61  V17.196  V18.110  V16.227   V14.47   V23.22    V2.12   V23.29 
  4.1761   3.8285   4.0208   3.9927   8.4033   9.2821   3.0228  14.6312 
  V12.43     R9.7   A157.3   V23.81   V23.82   V12.53   V23.83   V12.56 
 14.5555  17.0486  15.9813  19.0779  18.5811  18.7827  17.5376  20.5445 
 A152.84   V16.50  V22.122   V16.41    V4.32   V12.66  V19.245     V4.8 
 19.9792  19.8121  18.7319  22.9201  22.4598  20.7307  22.4718  22.0932 
 A180.15   V18.34  V20.213  V19.222  A180.39  V16.189   V12.18    V7.67 
 21.4827  23.3740  23.3249  22.8605  24.2542  25.6955  25.5307  23.2759 
 V17.165  V19.310  V16.190 A153.154  V19.308  V22.172   V10.98  V22.219 
 23.6758  23.0183  24.5248  25.3909  25.8264  26.3325  24.0536  25.4655 
  V16.33  V22.204  V20.167   V10.89   V12.79  V19.216   V14.90  A180.72 
 26.3604  25.8272  26.7432  26.4347  26.0909  25.7002  25.8305  26.3170 
  V16.21  A180.76  V15.164  A180.78    V14.5   V3.128  A179.13    V9.31 
 26.7791  26.6893  26.8161  25.9479  26.8964  26.7957  26.3988  26.0191 
 V20.230    V20.7  V20.234   V18.21  V12.122 
 26.6089  27.2027  26.7317  26.9673  26.8097 

Tuning?

With PCR, PLS, WA-PLS, and MAT we face an extra challenge

How to tune the model’s complexity parameter: \(l\) or \(k\)

How many components should we use? How many analogues?

Technically, we will underestimated RMSEP if we use CV to tune the model

We need a tuning test set

Tuning?

Form a test set from the training set (representative, evenly spaced) — set aside

Form a tuning or optimisation set from the training set (as above) — set aside

Build your TF using the remaining samples

1. Use CV to generate RMSEP for the tuning set sample
2. Find the value of \(l\) or \(k\) that gives lowest RMSEP on tuning set,
3. Use CV to generate RMSEP for the test set using the selected value of \(l\) or \(k\)

Close analogues

A measure of reliability for the reconstructed values can be determined from the distance between each fossil sample and the training set samples

For a reconstructed value to be viewed as more reliable, it should have at least one close modern analogue in the training set

Close modern analogues are defined as those modern training set samples that are as similar to a fossil sample as a low percentile of the observed distribution dissimilarities in the training set, say the 5\(^{th}\) percentile

Close analogues

Which samples in the core have a close match in the training set?

v12_mdc <- minDC(wa_m, V12.122 / 100, method = "chord")
plot(v12_mdc, use.labels = TRUE, xlab = "Depth")

Sample-specific errors

We can use the bootstrap approach to generate sample specific errors for each fossil sample

\[\mathrm{RMSEP} = \sqrt{s^2_{1_{fossil}} + s^2_{2_{model}}}\]

\(s^2_{1_{fossil}}\) is the standard deviation of the bootstrap estimates for the fossil samples

\(s^2_{2_{model}}\) is the average bias, the mean of the bootstrap OOB residuals from the model

Sample-specific errors

set.seed(1234)
v12_boot <- predict(wa_m, V12.122 / 100, bootstrap = TRUE, n.boot = 999)
reconPlot(v12_boot, use.labels = TRUE, ylab = "Summer SST", xlab = "Depth", display.error = "bars")

Sick science?

Juggins (Quaternary Science Reviews, 2013)

Are transfer functions too good to be true?

\(\mathbf{x}\) must have biological/physiological control on the taxa

Secondary nuisance variables are (very) problematic — they can bias the estimated species–environment relationships

You could be reconstructing variation in something else!

Variation downcore that has nothing to do with \(\mathbf{x}\) will show up as variation in \(\mathbf{x}_0\)

SWAP Example

Work through the SWAP example